const char* dgemm_desc = "Naive, three-loop dgemm.";

/*#include <memory.h> */
#include <stdlib.h>
#include <stdio.h>
#include <emmintrin.h>
/* This routine performs a dgemm operation
 *  C := C + A * B
 * where A, B, and C are lda-by-lda matrices stored in column-major format.
 * On exit, A and B maintain their input values. */    
void square_dgemm (int n, double* A, double* B, double* restrict C)
{
    /* Reorganize A to be row-major */
    /* TODO: How to do this efficiently? */
    double  *Atemp = (double *)malloc(sizeof(double)*n*n);
    for (int i = 0; i < n; ++i)
    {
        for (int j = 0; j < n; ++j)
        {
            Atemp[i*n+j] = A[i+n*j];
        }
    }
    A = Atemp;
  /* How to fix this so we actually free memory instead of seg faulting??? */
  /* free(Atemp);*/ 

// this only slowed things down
//   A = __builtin_assume_aligned(A,8);

   /*
   int x = __alignof__ (double);
   int y = __alignof__ (A);
   printf("Alignment of double is %d, Alignment of A is %d",x,y);
   */
   // Alignment of A and double are both 8!

  /* For each row i of A */
//  __m128d row,col,cij;
/* for (int i = 0; i < n; ++i)
    // For each column j of B 
    for (int j = 0; j < n; ++j) 
    {
      // Compute C(i,j) 
      __m128d cij = _mm_loadu_pd(C+i+j*n);
      double *a = A+i*n;
      double *b = B+j*n;
	for( int k = 0; k < n; k++ ){
	    __m128d row = _mm_load_pd(a+k);
	    __m128d col = _mm_load_pd(b+k);
	    cij = _mm_add_pd(cij,_mm_mul_pd(row,col));
	}
      _mm_storeu_pd(C+i+j*n,cij);
    }
}*/

__m128d c1 = _mm_loadu_pd(C+0*2);
__m128d c2 = _mm_loadu_pd(C+1*2);
for(int i = 0; i<2; i++){
a = _mm_load_pd(A+i*2);
b1 = _mm_load1_pd(B+i+0*2);
b2 = _mm_load1_pd(B+i+1*2);
c1 = _mm_add_pd(c1,_mm_mul_pd(a,b1));
c2 = _mm_add_pd(c2,_mm_mul_pd(a,b2));
}
_mm_storeu_pd(C+0*2,c1);
_mm_storeu_pd(C+1*2,c2);

free(Atemp);
